rmarkdown::render(‘./3_Clustering/3_Clustering.Rmd’)
Changes in myeloid and kidney cells after CLP - Analysis of 2 x 10X scRNA-seq samples from 2 pools of WT mice (3 Sham + 3 CLP): comparison of gene expression in different cell populations
indir <- "./processedData/2_1_Resolution_choice"
outdir <- "./processedData/3_Clustering"
dir.create(outdir, recursive = T)
library(Seurat)
integrated <- readRDS(paste0(indir, "/8.integrated.rds"))
integrated
## An object of class Seurat
## 24399 features across 18055 samples within 2 assays
## Active assay: integrated (2000 features, 2000 variable features)
## 1 other assay present: RNA
## 2 dimensional reductions calculated: pca, umap
Idents(integrated) <- "integrated_snn_res.1.7"
table(integrated@active.ident)
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
## 1310 1139 1102 1017 1006 914 831 811 763 715 698 642 528 504 456 446
## 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
## 431 425 417 352 348 346 334 327 238 207 204 192 184 144 142 133
## 32 33 34 35 36 37 38 39
## 128 101 92 90 88 88 81 81
pal <- colorRampPalette(c("#12999E", "#FAEB09", "#E82564", "#03539C"))
levels <- levels(integrated$integrated_snn_res.1.7)
colors.clusters <- pal(length(levels))
names(colors.clusters) <- levels
colors.clusters
## 0 1 2 3 4 5 6 7
## "#12999E" "#239F92" "#35A587" "#47AB7B" "#59B270" "#6BB864" "#7DBE59" "#8EC54D"
## 8 9 10 11 12 13 14 15
## "#A0CB42" "#B2D136" "#C4D82B" "#D6DE1F" "#E8E414" "#FAEB09" "#F8DB10" "#F7CC16"
## 16 17 18 19 20 21 22 23
## "#F5BD1E" "#F4AE25" "#F39E2C" "#F18F33" "#F08039" "#EE7141" "#ED6148" "#EC524F"
## 24 25 26 27 28 29 30 31
## "#EA4356" "#E9345D" "#E82564" "#D62868" "#C42C6C" "#B32F70" "#A13375" "#8F3679"
## 32 33 34 35 36 37 38 39
## "#7E3A7D" "#6C3D82" "#5B4186" "#49448A" "#37488F" "#264B93" "#144F97" "#03539C"
slices <- rep(1, length(levels))
pie(slices, col = colors.clusters, labels = names(colors.clusters))
d <- DimPlot(integrated, reduction = "umap", pt.size = 0.2, label = T,
label.size = 6, cols = colors.clusters)
d
pdf(paste0(outdir, "/1_DimPlot_umap_clusters_pc50_res1_7.pdf"),
width = 10, height = 8)
d
dev.off()
## png
## 2
colors.samples <- c("#12999E", "#FDA908")
names(colors.samples) <- levels(as.factor(integrated$sample.id))
slices <- rep(1, length(colors.samples))
pie(slices, col = colors.samples, labels = names(colors.samples))
p1 <- DimPlot(integrated, reduction = "umap", group.by = "sample.id",
pt.size = 0.2, cols = colors.samples)
p2 <- DimPlot(integrated, reduction = "umap", label = TRUE, pt.size = 0.2,
label.size = 6, cols = colors.clusters)
library(cowplot)
plot_grid(p1, p2)
pdf(paste0(outdir, "/2_2DimPlots_umap_samples_clusters_pc50_res1_7.pdf"),
width = 18, height = 8)
plot_grid(p1, p2)
dev.off()
## png
## 2
d <- DimPlot(integrated, reduction = "umap", group.by = "sample.id",
split.by = "sample.id", pt.size = 0.2, ncol = 2, cols = colors.samples)
d
pdf(paste0(outdir, "/3_DimPlot_umap_split_by_samples.pdf"), width = 16,
height = 9)
d
dev.off()
## png
## 2
f <- FeaturePlot(integrated, features = c("Nphs2", "Slc5a2",
"Clcnka", "Slc12a1", "Ptgs2", "Slc12a3", "Calb1", "Aqp2",
"Slc4a1", "Slc26a4", "Slc14a2", "Upk1a", "Cd22", "Adgre1",
"Pecam1", "Pdgfrb", "Cd68", "Cd14", "Acta2", "Csf3r", "Cd4"),
min.cutoff = "q9", order = T)
f
pdf(paste0(outdir, "/4_FeaturePlot_cellID.pdf"), width = 28,
height = 42)
f
dev.off()
## png
## 2
##Annotation of markers based on cluster markers from Susztak Science paper (Park et al., Science 360, 758–763 (2018) and Kidney International (2019) 95, 787–796; https://doi.org/10.1016/
https://science.sciencemag.org/content/360/6390/758.long https://www.kidney-international.org/article/S0085-2538(18)30912-8/fulltext
#Podocyte markers -> cluster 28
f2 <- FeaturePlot(integrated, features = c("Nphs2", "Podxl"),
min.cutoff = "q9")
f2
pdf(paste0(outdir, "/5_FeaturePlot_Podo.pdf"), width = 14, height = 7)
f2
dev.off()
## png
## 2
#Endothelial markers -> cluster 15
f3 <- FeaturePlot(integrated, features = c("Plat", "Pecam1"),
min.cutoff = "q9")
f3
pdf(paste0(outdir, "/6_FeaturePlot_Endo.pdf"), width = 14, height = 7)
f3
dev.off()
## png
## 2
#PT-S1 markers -> clusters 7,8,9
f4 <- FeaturePlot(integrated, features = c("Slc5a2", "Slc5a12"),
min.cutoff = "q9")
f4
pdf(paste0(outdir, "/7_FeaturePlot_PTs1.pdf"), width = 14, height = 7)
f4
dev.off()
## png
## 2
#PT-S2 markers
f5 <- FeaturePlot(integrated, features = c("Fxyd2", "Hrsp12"),
min.cutoff = "q9")
f5
pdf(paste0(outdir, "/8_FeaturePlot_PTs2.pdf"), width = 7, height = 7)
f5
dev.off()
## png
## 2
#PT-S3 markers -> cluster 5
f6 <- FeaturePlot(integrated, features = c("Atp11a", "Slc13a3"),
min.cutoff = "q9")
f6
pdf(paste0(outdir, "/9_FeaturePlot_PTs3.pdf"), width = 14, height = 7)
f6
dev.off()
## png
## 2
#Loop of Henle -> clusters 11, 13, 18
f7 <- FeaturePlot(integrated, features = c("Slc12a1", "Umod"),
min.cutoff = "q9")
f7
pdf(paste0(outdir, "/10_FeaturePlot_LOH.pdf"), width = 14, height = 7)
f7
dev.off()
## png
## 2
#Distal CT -> cluster 10
f8 <- FeaturePlot(integrated, features = c("Slc12a3", "Pvalb"),
min.cutoff = "q9")
f8
pdf(paste0(outdir, "/11_FeaturePlot_DCT.pdf"), width = 14, height = 7)
f8
dev.off()
## png
## 2
#Conn Tubule -> clusters 6, 20, 21, 29
f21 <- FeaturePlot(integrated, features = c("Calb1"), min.cutoff = "q9")
f21
pdf(paste0(outdir, "/12_FeaturePlot_ConnTub.pdf"), width = 7,
height = 7)
f21
dev.off()
## png
## 2
#CD PC -> cluster 21
f9 <- FeaturePlot(integrated, features = c("Aqp2", "Hsd11b2"),
min.cutoff = "q9")
f9
pdf(paste0(outdir, "/13_FeaturePlot_CD-PC.pdf"), width = 14,
height = 7)
f9
dev.off()
## png
## 2
#CD-IC -> clusters 24, 29, 39
f10 <- FeaturePlot(integrated, features = c("Atp6v1g3", "Atp6v0d2"),
min.cutoff = "q9")
f10
pdf(paste0(outdir, "/14_FeaturePlot_CD-IC.pdf"), width = 14,
height = 7)
f10
dev.off()
## png
## 2
#CD Trans -> cluster 29
f11 <- FeaturePlot(integrated, features = c("Slc26a4", "Insrr",
"Rhbg"), min.cutoff = "q9")
f11
pdf(paste0(outdir, "/15_FeaturePlot_CD-Trans.pdf"), width = 14,
height = 14)
f11
dev.off()
## png
## 2
#Fibroblast
f12 <- FeaturePlot(integrated, features = c("Plac8", "S100a4",
"Pdgfrb"), min.cutoff = "q9")
f12
pdf(paste0(outdir, "/16_FeaturePlot_Fib.pdf"), width = 14, height = 14)
f12
dev.off()
## png
## 2
#Macro -> cluster 22
f13 <- FeaturePlot(integrated, features = c("C1qa", "Cd68", "C1qb"),
min.cutoff = "q9")
f13
pdf(paste0(outdir, "/17_FeaturePlot_Macro.pdf"), width = 14,
height = 14)
f13
dev.off()
## png
## 2
#PMN -> cluster 36
f14 <- FeaturePlot(integrated, features = c("S100a8", "Ly6g",
"S100a9"), min.cutoff = "q9")
f14
pdf(paste0(outdir, "/18_FeaturePlot_PMN.pdf"), width = 14, height = 14)
f14
dev.off()
## png
## 2
#B lymph -> cluster 37
f15 <- FeaturePlot(integrated, features = c("Cd79a", "Cd79b",
"Cd19"), min.cutoff = "q9")
f15
pdf(paste0(outdir, "/19_FeaturePlot_Blymph.pdf"), width = 14,
height = 14)
f15
dev.off()
## png
## 2
#Tlymph -> cluster 30
f16 <- FeaturePlot(integrated, features = c("Ltb", "Cd4", "Cxcr6"),
min.cutoff = "q9")
f16
pdf(paste0(outdir, "/20_FeaturePlot_Tlymph.pdf"), width = 14,
height = 14)
f16
dev.off()
## png
## 2
#NK -> cluster 30
f17 <- FeaturePlot(integrated, features = c("Gzma", "Nkg7"),
min.cutoff = "q9")
f17
pdf(paste0(outdir, "/21_FeaturePlot_NK.pdf"), width = 14, height = 7)
f17
dev.off()
## png
## 2
#Novel1
f18 <- FeaturePlot(integrated, features = c("Slc27a2", "Lrp2",
"Cdca3"), min.cutoff = "q9")
f18
pdf(paste0(outdir, "/22_FeaturePlot_Novel1.pdf"), width = 14,
height = 14)
f18
dev.off()
## png
## 2
# library(Seurat)
DefaultAssay(integrated) <- "RNA"
clusters <- levels(integrated@active.ident)
conserved.markers <- data.frame(matrix(ncol = 14))
for (c in clusters) {
print(c)
markers.c <- FindConservedMarkers(integrated, ident.1 = c,
grouping.var = "sample.id", verbose = T)
markers.c <- cbind(data.frame(cluster = rep(c, dim(markers.c)[1]),
gene = rownames(markers.c)), markers.c)
write.table(markers.c, file = paste0(outdir, "/23_markers_",
c, ".txt"))
colnames(conserved.markers) <- colnames(markers.c)
conserved.markers <- rbind(conserved.markers, markers.c)
head(conserved.markers)
}
## [1] "0"
## [1] "1"
## [1] "2"
## [1] "3"
## [1] "4"
## [1] "5"
## [1] "6"
## [1] "7"
## [1] "8"
## [1] "9"
## [1] "10"
## [1] "11"
## [1] "12"
## [1] "13"
## [1] "14"
## [1] "15"
## [1] "16"
## [1] "17"
## [1] "18"
## [1] "19"
## [1] "20"
## [1] "21"
## [1] "22"
## [1] "23"
## [1] "24"
## [1] "25"
## [1] "26"
## [1] "27"
## [1] "28"
## [1] "29"
## [1] "30"
## [1] "31"
## [1] "32"
## [1] "33"
## [1] "34"
## [1] "35"
## [1] "36"
## [1] "37"
## [1] "38"
## [1] "39"
conserved.markers <- conserved.markers[-1, ]
openxlsx::write.xlsx(conserved.markers, file = paste0(outdir,
"/23_conserved.markers.xlsx"))
library(Seurat)
integrated <- RenameIdents(integrated, `5` = "PT-s3", `6` = "ConnTub",
`7` = "PT-s1", `8` = "PT-s1", `9` = "PT-s1", `10` = "DCT",
`11` = "LOH", `13` = "LOH", `15` = "Endo", `18` = "LOH",
`20` = "ConnTub", `22` = "Macro", `24` = "CD-IC", `28` = "Podo",
`36` = "PMN", `37` = "B lymph", `39` = "CD-IC")
d2 <- DimPlot(integrated, label = TRUE, label.size = 4)
d2
pdf(paste0(outdir, "/24_Dimplot_newidents.pdf"), width = 13,
height = 9)
d2
dev.off()
## png
## 2
d3 <- DimPlot(integrated, group.by = "sample.id", split.by = "sample.id",
pt.size = 0.2, ncol = 2)
d3
pdf(paste0(outdir, "/25_DimPlot_newidents_split_by_samples.pdf"),
width = 16, height = 9)
d3
dev.off()
## png
## 2
DefaultAssay(integrated) <- "RNA"
f19 <- FeaturePlot(integrated, features = "Il6", order = T, label = T,
label.size = 3)
f19
pdf(paste0(outdir, "/26_FeaturePlot_Il6.pdf"), width = 11, height = 10)
f19
dev.off()
## png
## 2
f20 <- FeaturePlot(integrated, features = c("Il6"), split.by = "sample.id",
max.cutoff = 3, cols = c("grey", "red"), order = T)
f20
pdf(paste0(outdir, "/27_FeaturePlot_Il6-sham-CLP.pdf"), width = 19,
height = 10)
f20
dev.off()
## png
## 2
library(ggplot2)
library(cowplot)
theme_set(theme_cowplot())
integrated$celltype.stim <- paste(Idents(integrated), integrated$sample.id,
sep = "_")
integrated$celltype <- Idents(integrated)
Idents(integrated) <- "celltype"
plots <- VlnPlot(integrated, features = c("Il6"), split.by = "sample.id",
group.by = "celltype", pt.size = 0, combine = FALSE)
library(patchwork)
wrap_plots(plots = plots, ncol = 1)
d <- DotPlot(integrated, features = "Il6", group.by = "celltype.stim")
openxlsx::write.xlsx(d$data, paste0(outdir, "/28_IL6_expn_per_celltype_stim.xlsx"))
d
pdf(paste0(outdir, "/29_DotPlot_IL6_celltype_stim.pdf"), width = 6,
height = 8)
d
dev.off()
## png
## 2
cluster33 <- WhichCells(integrated, idents = "33")
# others <- WhichCells(integrated, idents = "33", invert = T)
d <- DimPlot(integrated, label=T, group.by="celltype", cells.highlight= list(cluster33), cols.highlight = c("darkblue"
# , "darkred"
), cols= "grey")
d
pdf(paste0(outdir, "/30_DimPlot_integrated_label_group.by_celltype_cell.highlight_cluster33.pdf"))
d
dev.off()
## png
## 2
saveRDS(integrated, paste0(outdir, "/31.integrated.rds"))
sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-conda-linux-gnu (64-bit)
## Running under: Red Hat Enterprise Linux Server release 6.8 (Santiago)
##
## Matrix products: default
## BLAS: /gpfs/fs1/data/omicscore/Privratsky-Privratsky-20210215/scripts/conda/envs/privratsky/lib/libblas.so.3.8.0
## LAPACK: /gpfs/fs1/data/omicscore/Privratsky-Privratsky-20210215/scripts/conda/envs/privratsky/lib/liblapack.so.3.8.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] patchwork_1.1.1 cowplot_1.1.1 ggplot2_3.3.3 SeuratObject_4.0.0
## [5] Seurat_4.0.0
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-152 matrixStats_0.58.0 RcppAnnoy_0.0.18
## [4] RColorBrewer_1.1-2 httr_1.4.2 sctransform_0.3.2
## [7] tools_4.0.3 R6_2.5.0 irlba_2.3.3
## [10] rpart_4.1-15 KernSmooth_2.23-18 uwot_0.1.10
## [13] mgcv_1.8-33 lazyeval_0.2.2 colorspace_2.0-0
## [16] withr_2.4.1 tidyselect_1.1.0 gridExtra_2.3
## [19] compiler_4.0.3 formatR_1.7 plotly_4.9.3
## [22] labeling_0.4.2 scales_1.1.1 lmtest_0.9-38
## [25] spatstat.data_2.0-0 ggridges_0.5.3 pbapply_1.4-3
## [28] spatstat_1.64-1 goftest_1.2-2 stringr_1.4.0
## [31] digest_0.6.27 spatstat.utils_2.0-0 rmarkdown_2.6
## [34] pkgconfig_2.0.3 htmltools_0.5.1.1 parallelly_1.23.0
## [37] highr_0.8 fastmap_1.1.0 htmlwidgets_1.5.3
## [40] rlang_0.4.10 shiny_1.6.0 farver_2.0.3
## [43] generics_0.1.0 zoo_1.8-8 jsonlite_1.7.2
## [46] ica_1.0-2 zip_2.1.1 dplyr_1.0.4
## [49] magrittr_2.0.1 Matrix_1.3-2 Rcpp_1.0.6
## [52] munsell_0.5.0 abind_1.4-5 reticulate_1.18
## [55] lifecycle_1.0.0 stringi_1.5.3 yaml_2.2.1
## [58] MASS_7.3-53.1 Rtsne_0.15 plyr_1.8.6
## [61] grid_4.0.3 parallel_4.0.3 listenv_0.8.0
## [64] promises_1.2.0.1 ggrepel_0.9.1 crayon_1.4.1
## [67] deldir_0.2-9 miniUI_0.1.1.1 lattice_0.20-41
## [70] splines_4.0.3 tensor_1.5 knitr_1.31
## [73] pillar_1.4.7 igraph_1.2.6 future.apply_1.7.0
## [76] reshape2_1.4.4 codetools_0.2-18 leiden_0.3.7
## [79] glue_1.4.2 evaluate_0.14 data.table_1.13.6
## [82] vctrs_0.3.6 png_0.1-7 httpuv_1.5.5
## [85] gtable_0.3.0 RANN_2.6.1 purrr_0.3.4
## [88] polyclip_1.10-0 tidyr_1.1.2 scattermore_0.7
## [91] future_1.21.0 openxlsx_4.2.3 xfun_0.20
## [94] mime_0.10 xtable_1.8-4 later_1.1.0.1
## [97] survival_3.2-7 viridisLite_0.3.0 tibble_3.0.6
## [100] cluster_2.1.1 globals_0.14.0 fitdistrplus_1.1-3
## [103] ellipsis_0.3.1 ROCR_1.0-11
writeLines(capture.output(sessionInfo()), "./scripts/3_Clustering/3_Clustering.sessionInfo.txt")